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INTRODUCTION 


A  well  known  advantage  of  variational  solution  formulation  to 
boundary  value  problems  la  that  the  differentiability  requirements  of  the 
approximate  solutions  can  be  relaxed.  For  initial  value  problems,  however, 
this  advantage  is  somewhat  diminished  by  the  complication  due  to  the 
appearance  of  the  farther  end  condition.  This  complication  can  be  eliminated 
by  the  use  of  an  adjoint  variational  principle  as  we  have  demonstrated  for  a 
simple  initial  value  problem  in  a  previous  paper. 1  The  more  involved  anal¬ 
ysis  for  the  mixed  initial-boundary  value  problem  has  also  been  worked  out. 

The  purpose  of  this  report  is  to  employ  the  adjoint  variational  principle 
in  the  form  of  finite  element  formulation  for  solving  the  stress  wave 
problems.  The  hyperbolic  partial  differential  equation  governing  the  motion 
is  second  order  both  in  spatial  and  time  domains. 

Ly(x,t)  +  Q(x,t)  -  0  (1) 

where 

Ly  -  (ayt)t  +  (*yx)x  (2) 

Ue  seek  explicitly  the  numerical  transient  solutions  of  y,  y^,  yx  and  yxt  for 
assigned  boundary  and  initial  conditions.  The  term  yx  will  give  the  stress 
wave  in  a  longitudinal  bar.  The  study  is  the  extension  of  previous  work  on 
initial  and  boundary  value  problems. 1 


1-Shen,  C.  N. ,  "Method  of  Solution  For  Variational  Principle  Using  Bicubic 
Hermlte  Polynomial,"  presented  at  the  27th  Conference  of  Army 
Mathematicians,  West  Point,  NY,  June  1981. 


INTEGRAL  OF  BILINEAR  EXPRESSION 

The  integral  of  a  bilinear  expression  for  a  two-dimensional  problem 
having  second  order  partial  derivatives  in  both  space  and  time  can  be  written 

as 

Xb  fcb 

I  ■  /  /  »i[y(x  ,t)  ,y(x»t)  Jdtdx  (3) 

Xo  to 

where  ft[y,y]  is  a  given  bilinear  expression  in  the  form 

n[y,y]  -  <*ytyt  +  *yxyx  (*) 

The  quantity  y  is  the  adjoint  of  y  and  the  subscripts  t  and  x  Indicate  the 
partial  derivatives  of  the  functions  y  and  y. 

Two  different  forms  of  integrals  and  end  conditions  can  be  obtained  from 
Eq.  (A).  The  first  form  is  obtained  by  integrating  by  parts  on  the  adjoint 
variable. 

xb  tb  -  xb  -  tb  tb  Xb 

I  ■  -/  /  yLydtdx  +  /  ayty|  dx  +  /  iyxy|  dt  (5) 

Xo  to  Xo  to  to  Xo 

where  Ly  is  given  in  Eq.  (2).  ✓ 

In  addition,  we  can  perform  Integration  on  the  original  variable  to  give 

xb  tb  —  xb  -  tb  tb  xb 

I  -  -/  /  yLydtdx  +  /  ayty|  dx  +  /  4yxy|  dt  (6) 

Xo  to  *o  to  to  Xq 

where 

Ly  -  (ayt)t  +  (*yx>x  (7) 

In  a  previous  paper *  we  show  that  the  bilinear  concomitant  D  has  to  be 
identically  zero,  l.e.. 


Ishen,  C.  N. ,  "Method  of  Solution  For  Variational  Principle  Using  Bicubic 
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Mathematicians,  West  Point,  NT,  June  1981. 
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D  ■  /  /  yLydCdx  -  /  /  yLydtdx 

*o  tQ  *0 


(8) 


By  equating  Eqs.  (5)  and  (6)  and  solving  for  D  in  Eq.  (8),  we  are  converting 
the  double  Integral  into  two  single  integrals  in  terms  of  the  initial  and 
boundary  conditions. 

We  can  express  the  quantity  D  as  the  sum  of  two  parts  for  end  conditions 
as  Dj  and  D2«  Thus  one  defines 

D  -  Dj  +  D2 

The  terms  in  Dj  Involve  the  initial  conditions  of  y  and  y  as 


xb  -  tb  -  tb 
D1  •  /  Wtyl  -  ayty|  )dx 

xO  t0  tQ 


*b 


■  /  ^b(ytbyb-ytbyb)  -  “o(ytoyo-ytoyo>}<ix 

*o 

The  terms  in  D2  involve  the  boundary  conditions  of  y  and  y  as 


tb  -  xb  -  xb 
02  -  J  Uyxy|  -  lyxy|  >dt 

tQ  Xq  Xq 


^b 


I  t *b ( yxbyb-yxbyb )  -  *o(yxoyo-yxoyo)><it 

to 


In  order  that  D  =  0  in  Eq.  (9)  it  is  sufficient  that 

D!  =  0 
D2  h  0 


and 


(9) 


(10) 


(ID 


(12a) 

(12b) 


END  CONDITIONS  FOR  THE  ADJOINT  SYSTEMS 

In  oruiir  to  satisfy  the  two  requirements  in  Eq.  (12),  we  separate  them  in 


two  parts.  Let  us  consider  first  the  time  domain  and  assume  that  the  adjoint 

variables  are  assigned  as 


o 

>> 

i 

£ 

»  yo  -  yb 

ytb  -  -<»b"lao  yto 

.  yto  -  ytb 

jr; 

ab  *  0 

<*o  *  0 

The  above 

adjoint 

initial  conditions 

satisfy  the  requirement 

tha  t  D , 

Eq.  (10). 

Now  we 

turn  to  the  spatial 

domain  and  assume  that 

the  ad 

-  ? 

variables 

are 

yb  -  6yb 

Yo  •  8y0 

yxb  -  Syxb 

Yxo  "  ^Yxo 

r1 

The  above  adjoint  boundary  conditions  satisfy  the  requirement  that 
Eq.  (11),  with  0  as  a  constant. 

By  giving  the  appropriate  values  of  these  adjoint  variables  in 

the  original  variables,  one  nay  find  that  the  requirement  D  =  0  can 

satisfied.  This  leads  to  the  result1-  previously  found  as 

*b  *b  -  tb  *b  - 

Jly.y]  mJ  J  Qydtdx  +  J  /  y(QH.  dtdx  -  o 

to  *0  to  Xo 

FIRST  VARIATION 

By  taking  variation  on  Eq.  (18)  we  have 

6J  -  <5J[6y]  +  6 J[ 6y ] 


tb  Xb  —  tb  Xb  - 

-  /  /  6y(Ly)dtdx  +  /  /  y(L6y)dtdx 

to  Xq  to  Xq 


-  0 
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6j[6y]  ■  J  /  dy(Ly+Q)dtdx 

t0  *o 

tb  *b - 

«Jl«y]  -  /  /  6y(Ly+Q)dtdx 

t0  Xq 


Since  D  s  0  In  Eq.  (8)  the  variation  6D  should  be  zero 

50  ■  6D[6y]  +  6D[5y]  -  0  (22) 

Since  the  variations  6y  and  6y  are  independent,  then 

tb  Xb  -  tb  Xb 

6 D[ 6y ]  -  /  J  y(L6y)dtdx  -  J  J  5y(Ly)dtdx  -  0  (23) 

to  Xo  to  Xq 

which  la  the  same  as  the  last  two  terms  In  Eq.  (19)  which  vanish.  Thus 

5 J  -  6J[6y]  +  6J[6y]  -  0  %  (24) 

Since  the  variations  5y  and  5y  are  independent 

tb  «b  * 

6j[5y]  -  /  /  6y(Ly+Q)dtdx  -  0  (25) 

to  xo 

where  Ly  is  given  in  Eq.  (2)  and  contains  higher  derivatives  than  the  first 
partlals  in  y.  It  is  Intended  to  Include  only  lower  order  partial 
differentiation  in  y.  This  can  be  achieved  by  considering  the  variations  of 
the  bilinear,  expression  I  given  by  Eqs.  (3)  and  (4)  as 

5l[Sy]  -  /  /  [aytSyt  +  AyxSyxMtdx  (26) 

to  *o 

A  different  form  of  the  above  variation  can  be  obtained  from  Eq.  (5)  as 


*b  tb  -  *b  -  tb  tb  -  xb 

6l[5y]  -  -/  /  SyLydtdx  +  /  5yoyt|  dx  +  /  dyiyxl  dt  (27) 


Equating  Eqs.  (26)  and  (27),  solving  fo  the  term  containing  integral  for  6yLy 

and  substituting  into  Eq.  (25)  we  have 

xb  -,tb  *b  -xb 

6J[6y]  -  /  ayt5yj  dx  +  j  Ayx «Sy |  dt 
x0  t0  to  xo 

Xb  tb  -  Xb  tb 

+  /  /  SyQdtdx  -  /  /  [ayt6yt  +  iyx6yx]dtdx  -  0  (28) 

Xq  tQ  Xq  t0 

This  is  the  key  equation  which  uses  variational  principle  in  solving  a  mixed 
initial  and  boundary  value  problem  for  a  wave  equation. 

DISCUSSION  OF  THE  VARIATIONAL  EQUATION 

Let  us  discuss  the  various  terms  in  Eq.  (28),  the  variational  formulation 
of  the  wave  equation,  into  three  parts  as  follows. 

(1)  The  initial  conditions  of  the  original  variables  are  given  and 
variations  of  the  adjoints  at  the  far  end  are  zero.  The  first  term  in  Eq. 

(28)  contains  the  product  of  yt^y  evaluated  at  the  initial  condition  yto^yo 
and  at  the  final  condition  ytb^yb*  Since  the  value  of  yb  is  known  by  Eqs. 

(13)  and  (16),  6yb  -  0.  That  is,  the  variations  of  the  adjoint  variable  at 
the  far  end  are  zero. 

(2)  The  boundary  conditions  of  the  original  variables  and  variation  of 
the  adjoints  can  be  determined.  The  second  term  in  Eq.  (28)  is  the  boundary 
term  involving  the  variation  5y  and  the  variable  yx.  For  a  longitudinal  or  a 
torsional  bar  the  end  conditions  are 

from  Eq.  (16)  for  the  fixed  end 


from  Eq.  (17)  for  the  free  end 

yx  ■  0  yx  •  0  6yx  ■  0  (30) 

The  variations  in  the  adjoint  variables  shown  in  the  last  column  coincide 
with  the  same  end  conditions  in  the  original  variables  given  in  the  first 
column,  whether  they  are  on  the  left  or  the  right  boundary. 

(3)  Interior  Region  -  The  last  two  terms  in  Eq.  (28)  give  the  interior 
region  where  the  forcing  function  Q,  the  adjoint  variations  6y,  6yt,  and  6yx 
and  the  variables  yt,  and  yx  are  shown.  No  second  order  partial  of  y  with 
respect  to  x  is  present.  Thus  the  variables  that  are  needed  for  the 
computation  are  y,  yt,  and  yx.  This  requires  a  c*  continuity  in  both  spatial 
and  time  domain. 

TRANSFORMATION  OF  COORDINATES 


The  Integral  signs  in  Eq.  (28)  can  be  converted  into  summation  signs  if 
discrete  Intervals  for  integration  are  used.  We  may  take  some  scale  factors 
to  nondlmensionalize  the  problem  by  giving 


tQ  “  0 

»  tb  “  1 

0  <  t  <  1 

(31) 

*o  -  0 

.  *b  ”  l 

0  <  x  <  1 

(32) 

Moreover,  Eq.  (28)  can  be  discretized  by  letting 

5  -  Ht-i+1 

0  <  5  <  1 

1  -  1,2, ... ,H 

(33) 

n  -  Kx-j+1 

0  <  n  <  l 

i  -  1,2 . K 

(34) 

where  H  and  K  are  number  of  intervals  for  t  and  x  respectively.  Thus  the 
partial  derivatives  are: 


-  ^ 

L'- 

r 


r 


r 


r  . 


7 


(36) 


3y 


3y 

an 


Kyn 


Use  of  Sqs.  (28),  (31)  through  (36)  then  leads  to 

0  -  <5 J[6y] 

K  h  1  cb 

m  l  -  f  ayr(^»j )6y(i»j )dn| 
j-1  K  0  tQ 

H  V  1  XK 

+  l  -  /  Ayn(1*j)6y(i»j)d5| 

1-1  H  0  x0 

K  H  i  11- 
+  I  I  —  /  /  Q6y(i» J )d^dn 

j-1  t-i  HK  0  0 

“I  1  t”  /  /  ay^(i»j)6y5(i»j)dCdn  +  -  /  /  lyn(i » j )6yn(* * j )d£dn} 
j-1  1-1  K  0  0  H  0  0 

(37) 

SPLINE  FUNCTION 

We  may  express  the  variables  y(^*J)  and  <5y(i»j)  in  Eq.  (37)  In  terms  of 
the  (1x16)  spline  function  aT(5,n)  and  the  (16x1)  node  point  function  y(*»J) 
as  follows. 

y(i»J)(£,n)  -  aT(C,h)Y(i»3)  (38) 

where 

aT(£,u)  -  { [a1 (^ *n) ]T  la2(C,n)!T  [a3(S,n)]T  [a4(E,n)lT}  (39) 

and 

6y(1*J)(t,n)  -  aT(t,n)6Y(i.3)  (40) 


A  typical  term  for  a  product  can  be  written  as 

-  [6Y(1»J ) ]^a( £ ,n)a^(5 ,n)Y(l» J )  (41 

Thus  Eq.  (37)  becomes 

K  -  T 

6J(6y)  -  I  6Y(tb.j)]  Po5(tb)Y(tb.J) 
j-1 

K  ij 

-  I  «YCt0,j)  ]  P05(to)Y(to.j) 

J-1 

H  _  X 

+  l  [6Y<i.*b)]  Pon(xb)Y(i.*b> 
i-1 

H  T 

-  I  6Y<l.*o>l  P0n(xo)Y(i.*o) 
i-1 

K  H  X 

+  l  l  [6YU.J>]  q(i.J) 

J-1  i-1 

-  I  l  [6Y(1»J)]TP<1»J)Y(i*J)  -  0  (42 

J-1  i-1 

where  the  coefficient  P  contains  integrals  involving  the  spline  functions 
a(£,ri)  and  its  partial  derivatives  as  given  in  a  previous  paper. 1 


GRID  SYSTEMS  FOR  FINITE  ELEMENT 

We  take  a  finite  element  represented  by  the  (16x1)  vector  y(*»J)  which 
has  a  grid  of  four  (4x1)  vectors  Yj^^J)  through  Y^lfJ),  thus 

YU.J)  -  ( IY1C±»J> JT  [Y2(i,J))T  [y3(1»J))T  [Y4C1-*J  )]T>  (43; 
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Each  of  the  (4x1)  vectors  has  four  components,  consisting  of  the  function,  its 
first  partlals  in  both  directions,  and  its  mixed  partial,  as  shown  in  Figure 


1. 


These  vectors  are, 


yUi.n-j) 

yUi.nj+i) 

yc(5i»nj) 

Y3(i»j)  - 

ycCSi.hj+i) 

yn^i.^j) 

yn(£i>hj+i) 

ysn^i.hj) 

y5n(Si*nj+i> 

y(5i+l.nj) 

y(^i+l»nj+i) 

y£(Ci+l.nj) 

Y4U»j>  - 

yt^i+i.nj+i) 

yn(Si+i»hj) 

ynC^i+l.^j+l) 

ycn<^i+l.r'j) 

y^nC^i+l.^j+l) 

We  use  the  vertical  direction  for  the  temporal  domain.  If  we  increase  the  row 
index  from  i  to  i+1 ,  then  the  grid  point  shifts  down  by  one  step  and  the 
following  holds 

Yl(i+l,j)  .  Y2<1*J)  Y3(i+1tj)  -  Y^M)  (45) 

If  we  increase  the  column  index  from  j  to  j+1,  then  the  grid  point  shifts  to 
the  right  by  one  step  and  one  obtains 

YiUJ+l)  -  Y3^»J )  Y2(i»3+1>  -  Y4C 1 » J  >  (46) 

Figure  2  shows  the  relationship  of  the  grid  system  by  assembly  of  finite 
elements  in  the  horizontal  direction,  which  is  in  the  spatial  domain. 


ASSEMBLY  OF  MATRICES 

In  order  to  solve  Eq.  (42)  by  finite  element  method,  assembly  of  matrices 
from  local  form  into  global  form  is  necessary.  For  Instance,  the  last  term  of 
Eq.  (42)  is  taken  as  -6Jp(6y).  Then 

^  K  H 

6Jp(6y)  -  l  l  16Y<i.j)]Tp(1»^Y(i»j)  (47) 

j-1  i-1 

Since  we  know  that  the  internal  in  time  can  be  made  as  small  as  possible,  with 
H  -  1,  we  have 

K 

6Jp(6y)  -  I  [  6Y<  1  )  ]Tp(  1  )yO  > 

j-I  _  _  __ 

K 

-  I  {[6Yi(1,j)]TI6Y2(l,J)]T[6Y3(l.j)]T[6y4(l,J)]T}  pu  p12  p13  pu  Yl<l 
j-1 

**21  p22  p23  P24  Y2^»J) 


O.j) 


P31  p32  p33  p34  *3 


P41  p42  p43  p44 1  Y4(1»J> 


It  is  noted  from  Figure  2  that  the  variables  can  be  indexed  as 


Y3d»J)  -  YjOtj+l)  -  Y2j+i 
Y4<1»J)  -  Y2<l*J+D  -  Y2j+2 


j-0,l,...k 


For  j  -  0, 


For  j  -  k  -  5 


Yi<M>  -  Yi  ,  Y^1*1)  -  Y2 

Y3(1|S)  m  Yii  ,  Ya^*^)  -  Yj2 


(53) 


Also  from  Figure  2,  Che  adjoint  variations  are 


-  6Y1(1»J+1>  -  <$Y2j+i 

6Y4^1»J)  -  5Y2<l*j+1>  -  6Y2J+1 

(53) 

(54) 

For  j  -  0, 

dY^1*!)  -  6Yj  ,  6Y2(1»1)  -  6Y2 

(55) 

For  j  -  k  -  5, 

fiY3(1.5)  .  6Yn  ,  6Y4<1*5>  -  6Y12 

(56) 

Now  the  local  matrices  In  Eq.  (48)  can  be  assembled  into  a  global  band  matrix 
shown  In  Figure  3.  Those  elements  not  explicitly  written  are  zeroes  in  Figure 
3. 


Since  the  adjoint  variable  yj)  at  the  far  end  is  assigned  in  terms  of  the 
known  Initial  value  y0>  the  variation  is 

fiyb  -  0  (57) 

From  Figure  2  we  have 

6Y2  •  6Y4  -  6Yg  -  6Y8  ■  6Y10  ■  $Ti2  *  fi*EVEN  “  0  (58) 

This  Is  equivalent  to  deleting  the  even  rows  of  the  matrix  in  Figure  3.  The 
deletion  Is  marked  in  Figure  4.  The  number  of  relationships  is  reduced  to 
half  of  the  original  dimension. 

The  variables  Yodq  In  Figure  2  are  the  initial  values  of  the  problem 
which  are  supposed  to  be  given.  Thus,  Yj,  Y3,  Y5,  Y7,  Y9,  Yu  are  all  knowns. 
The  coefficients  related  to  these  knowns  should  eventually  be  shifted  to  the 


right  side  of  the  equation 


FURTHER  DELETIONS  AND  KNOWNS 


Suppose  we  have  a  bar  with  the  fixed  end  at  the  left.  Then  from  Eq.  (29) 


one  obtains 


o 

■ 

£ 

(59) 

and 

_ 

6y0  -  0 

(60) 

The 

above  equations 

translate  to  be 

y(2,l)  m  0  (known) 

(61) 

and 

6y(l|l)  »  0  (deletion) 

(62) 

On  i 

the  other  hand  we 

!  have  a  free  end  at  the  right.  Then  Eq. 

(30)  gives 

yXb  ■  o 

(63) 

and 

mm 

fiyxb  -  o 

(64) 

The 

above  equations 

yield 

' 

yn<2.6)  -  o  (known) 

(65) 

and 

6yn(l»6)  ■  0  (deletion) 

(66) 

Figure  S  gives  the  variation  of  adjoint  variables.  It  shows  two  extra 
zero  variations  at  the  first  row,  dy(l>D  at  the  left  and  Gy^1*6^  »t  the 
right.  We  have  also  all  zero  variations  on  the  second  row.  Figure  6  shows 
the  known  and  unknown  variables.  There  are  two  extra  known  variables  In  the 
second  row  due  to  boundary  conditions,  y(2»l)  at  the  left  and  yy/^»^  at  the 
right.  The  first  row  gives  all  known  Initial  conditions. 

CONCLUSIONS 

Direct  computation  of  stress,  l.e.,  numerical  solution  for  first  spatial 
derivatives  of  the  displacement  can  be  obtained  directly.  This  is  Important 
if  the  problem  has  noisy  components  in  the  solution  of  the  displacement. 


Computation  can  ba  made  successively,  i.e.,  the  final  values  of  the  solution 
at  the  first  stage  in  time  can  be  used  as  the  initial  values  of  the  second 
period  in  time.  The  variations  of  the  adjoint  variables  at  the  far  end  in 
time  for  an  Initial -boundary  value  problem  are  zeroes.  Deletion  of  many  rows 
in  the  assembled  matrix  is  possible.  The  assembled  matrix  for  computing  is 
reduced  to  less  than  half  size  in  linear  dimension,  from  (2nx2n)  to  (nxn). 
Hence,  a  bigger  number  of  Intervals  in  the  spatial  dimension  can  be  handled. 
The  reduced  matrix  is  a  band  matrix  which  makes  the  storage  requirement  for 
computation  much  easier. 
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Figure  3.  Matrix  Assembly  Global  Form 
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Figure  5.  Zero  Variation  of  Adjoint  Variables. 


Figure  6.  Variables,  Known  or  Unknown. 
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